
\documentclass[12pt]{article}
\usepackage{geometry} 
\usepackage{graphicx}
\usepackage{subfig}
\usepackage{bbold}
\geometry{a4paper} 
\title{OpenKet Brief Manual}
\author{V\'{i}ctor S C Canela}
\date{}

\begin{document}

\maketitle
\newpage
\tableofcontents
\newpage
\section{Introduction}
	This is an introductory manual to the Python library OpenKet. This library is 
	a tool for manipulating objects from quantum mechanics, such as vectors in the form 
	of kets, operators, etc. 

	It is free software under GNU-GPL v3 license, and uses other
        free software libraries, such as Sympy, Numpy and Scipy.

	The project started at the end of 2009 with Vicente
        Rodr\'{i}guez under the supervision of Pablo
        Barberis-Blostein\footnote{at IIMAS, UNAM}, as part of the
        requirements for his graduation.  I started working with Pablo
        one year later, and continued programming OpenKet.

	We want OpenKet to be as transparent as possible. This means
        making the input process close to writing it on paper. Dirac's
        notation is used to represent vectors as kets (having the same
        visual representation $|x\rangle$) for example.
	
	The basic structure of the library is the definition of object
        classes with specific rules regarding the basic algebra they
        should follow. Then there are functions, intended for a
        specific purpose, for example generate time evolution
        equations from a Hamiltonian.

	OpenKet can be downloaded \newline from
        http://code.google.com/p/openket/. In order to run OpenKet it
        is needed: Sympy, Pylab and Scipy\footnote{Working in iPython
          with Pylab is the ideal environment for OpenKet}).
	
        OpenKet is in its early stages and has a lot of bugs. Consequently
        is continuously improving. 


	Before typing the examples presented in the following
        sections, it is necessary to import OpenKet\footnote{if
          working in iPython type \texttt{run -i openket.py}} \texttt{
	\begin{tabbing}
	>>> from openket import *
	\end{tabbing}}

\newpage
\section{Kets and Bras}

	\subsection{Properties}
		In quantum mechanics, states are represented as vectors of a Hilbert 
		space. Using Dirac's notation, these vectors are represented as $|x\rangle$ and
		are called \emph{kets}. OpenKet understands kets as an object (class)
		with specific properties or methods. We can define a ket as follows, 
		
		\texttt{
		\begin{tabbing}
			>>> x = var("x") \\
			>>> v = Ket(x); v \\
			|x> \\
		\end{tabbing}}
		
		In this example \texttt{x} and \texttt{v} are two different objects. \texttt{x} 
		is a \emph{Sympy} variable and \texttt{v} is the object Ket. The argument given 
		for this object is the vector tag.
		This tag must be either a number class or 
		a \emph{Sympy} variable, as in the example. \\
		
		There is another argument (optional) which further specifies the ket. This is the 
		operator tag. The ket is then the eigenvector of the operator (whose name is 'operator tag'), 
		with eigenvalue equal to the 'vector tag'.
		
		\texttt{
		\begin{tabbing}
			>>> u = Ket(1,"Op1") \\
			>>> w = Ket(1, "Op2") \\
		\end{tabbing}}
		 
		Although \texttt{u} and \texttt{w} have the same vector tag, they are not treated 
		as the same object. For this kind of tags, strings as well as Sympy variables 
		may be used. \\
		
		The object Ket has 2 methods: \texttt{.op}\quad and \texttt{.eig}\quad used to extract all 
		the information from it. While \texttt{.eig}\quad returns the vector tag, \texttt{.op}\quad 
		returns the operator tag as follows,
		
		\texttt{
		\begin{tabbing}
			>>> u.eig \\
			1 \\
			>>> u.op \\
			'Op1'
		\end{tabbing}}
		
		Vectors corresponding to the dual space are represented by $\langle y|$ and are 
		called \emph{bras}. There is also an object Bra in OpenKet, with exactly 
		the same characteristics as de Ket object,
		
		\texttt{
		\begin{tabbing}
			>>> y = var("y") \\
			>>> z = Bra(y); z; z.eig; z.op \\
			<y| \\
			y \\
			'default' \\
		\end{tabbing}}

		Note that when no operator tag is specified \texttt{.op} returns \texttt{'default'}.
	
	\subsection{Basic Algebra}
		Basic operations may be performed on the Kets and Bras objects, such as:
		
		\texttt{
		\begin{tabbing}
			>>> w = Ket(0); x = Ket(1); y = Bra(0); z = Bra(1) \\
			>>> (2.1 + 2)*w + (4.3 + I)*x \\
			4.1|0> + (4.3 + I)|1> \\
			>>> y - 1.7*y + 0.5*z + I*z \\
			- 0.7<0| + (0.5 + I)<1| \\
		\end{tabbing}}
		
		Where the \texttt{I} $\equiv i = \sqrt{-1}$. Multiplication is also supported for example, 
		if we multiply a bra by a ket we get a number,
		
		\texttt{
		\begin{tabbing}
			>>> y*w \\
			1 \\
		\end{tabbing}}
		
		OpenKet assumes \texttt{w} and \texttt{y} are normalized. However, if we were to 
		define state vectors with norm not equal to unity, the results would be different,
		
		\texttt{
		\begin{tabbing}
			>>> psi = 2.1*Ket(0) + I*Ket(1); phi = I*Bra(0) - 4.5*Bra(1) \\
			>>> phi*psi \\
			-2.4*I \\
		\end{tabbing}}
		
		Multiplying a ket by a bra results in an outer product,
		
		\texttt{
		\begin{tabbing}
			>>> w*y; psi*phi \\
			|0><0| \\
			2.1*I|0><0| - 9.45|0><1| - |1><0| - 4.5*I|1><1| \\
		\end{tabbing}}
		
		Tensor products can be simulated having two different operator tags. For example,
		
		\texttt{
		\begin{tabbing}
			>>> f = Ket(0, "Op1"); g = Ket(0, "Op2"); h = Ket(1, "Op3") \\
			>>> l = Bra(1, "Op1"); m = Bra(1, "Op2"); n = Bra(1, "Op3")\\
			>>> f*g*h\\
			|0\_Op1>|0\_Op2>|1\_Op3>\\
			>>> l*m*n\\
			<1\_Op3|<1\_Op2|<1\_Op1|\\
		\end{tabbing}}

\section{Operators}
	\subsection{Properties}
		An operator in OpenKet is an object (class) defined as follows,
		
		\texttt{
		\begin{tabbing}
			>>> A = Operator("A"); A \\
			A \\
		\end{tabbing}}
		
		The symbol on the left of the equality does not need to be the same as the 
		argument of 
		\texttt{Operator}, this is just how you are calling it, whereas the argument is how 
		OpenKet tags it.\\
		 
		This can be made clear with the only method of \texttt{Operator}, \texttt{.op}. This 
		has the same effect as with kets and bras, it returns the tag of the operator 
		for example,
		
		\texttt{
		\begin{tabbing}
			>>> name = Operator("up"); name.op \\
			'up' \\
		\end{tabbing}}
		
		It is convenient to have the same tags if no confusion is produced. 
		However the option is always available to name objects the way we want.
		
	\subsection{Creation / Annihilation}
		In many physical systems it is convenient to make a canonical transformation from 
		the position - momentum operators to the creation - annihilation operators. The 
		obvious example is the quantum harmonic oscillator. In terms of these particular 
		operators, the Hamiltonian may be written,
		
		\[
		\hat{\mathcal{H}} = \hbar\omega\left(\hat{a}^{\dagger}\hat{a} + \frac{1}{2}\right)
		\]
		
		The number operator is defined as $\hat{N} = \hat{a}^{\dagger}\hat{a}$ so the above
		expression becomes,
		
		\[
		\hat{\mathcal{H}} = \hbar\omega\left( \hat{N} + \frac{1}{2}\right)
		\]
		
		The eigenvectors of $\hat{N}$ are called Fock states and therefore satisfy the equation,
		
		\[
		\hat{a}^{\dagger}\hat{a}|n\rangle = n|n\rangle
		\]
		
		Where $n$ is the energy level.\\
		Individually, the operators $\hat{a}^{\dagger}$ and $\hat{a}$ operate on $|n\rangle$ 
		as follows,
		
		\[
			\hat{a}^{\dagger}|n\rangle = \sqrt{n+1}|n+1\rangle \\ 
		\]		
		\[
			\hat{a}|n\rangle = \sqrt{n}|n-1\rangle
		\]
		
		OpenKet has two separate object classes for these operators. They are the 
		\texttt{CreationOperator} and the \texttt{AnnihilationOperator}. Both have 2 
		arguments (optional). The first is an 'operator tag' type; necessary when having
		tensor products of vectors with tags that could be confused. The second is a maximum 
		value of the ket, i.e. a number such that if the ket has that eigenvalue, applying
		a creation operator on it results in zero. This all looks like this,
		
		\texttt{
		\begin{tabbing}
			>>> a = AnnihilationOperator(); aa = CreationOperator(); n=var("n") \\
			>>> s = Ket(n); t = Bra(n) \\
			>>> a*s; aa*s; t*a; t*aa \\
			n**(1/2)|-1 + n> \\
			(1 + n)**(1/2)|1 + n> \\
			(1 + n)**(1/2)<1 + n| \\
			n**(1/2)<-1 + n| \\
		\end{tabbing}}
		
		More complex,
		
		\texttt{
		\begin{tabbing}
			>>> a1 = AnnihilationOperator("mode1") \\
			>>> aa1 = CreationOperator("mode1", 4) 
			\textsl{\# maximum value = 4} \\
			>>> a2 = AnnihilationOperator("mode2") \\ 
			>>> aa2 = CreationOperator("mode2", 9) 
			\textsl{\# maximum value = 9} \\
			>>> p = Ket(2, "mode1")*Ket(3, "mode2"); p \\
			|2\_mode1>|3\_mode2> 
			\textsl{\# bipartite system} \\
			>>> a1*p; a2*p; a1*a2*p \\
			2**(1/2)|1\_mode1>|3\_mode2> \\
			3**(1/2)|2\_mode1>|2\_mode2> \\
			6**(1/2)|1\_mode1>|2\_mode2> \\
			>>> aa1*a1*p; aa2*a2*p \\
			2|2\_mode1>|3\_mode2> 
			\textsl{\# number operator acting on 'mode1'}\\
			3|2\_mode1>|3\_mode2> 
			\textsl{\# number operator acting on 'mode2'}\\
		\end{tabbing}}
		
		It is possible to apply the operators as many times as we like. For example 
		(maintaining names),
		
		\texttt{
		\begin{tabbing}
			>>> (aa**4)*s \\
			(1 + n)**(1/2)*(2 + n)**(1/2)*(3 + n)**(1/2)*(4 + n)**(1/2)|4 + n> \\
			>>> (aa1**2)*p; (aa1**3)*p \\
			2*3**(1/2)|4\_mode1>|3\_mode2> \\
			0 \\
		\end{tabbing}}
		
		The last zero being a result of having surpassed the eigenvalue 4 (maximum value) for the 
		\texttt{photon1} system.
		
	\subsection{Examples}
		\subsubsection{Eigenvalue Syntaxis}
			The eigenvalue equation,
			
			\[
			\hat{A}|a\rangle = a|a\rangle
			\]
			
			Can be computed using OpenKet. It is necessary to match the names of the operator 
			tags of both the \texttt{Ket} and the \texttt{Operator},
			
			\texttt{
			\begin{tabbing}
				>>> a = var("a"); A = Operator("A"); u = Ket(a, "A") \\
				>>> A*u \\
				a|a\_A> \\
			\end{tabbing}}
			
		\subsubsection{Quantum Gates}\label{quantumgates}
			In quantum computing the basic unit of information is called a qubit, and is represented 
			as a linear combination of two vectors: $\{|0\rangle, |1\rangle\}$. 
			(computational basis). Manipulation of qubits 
			is possible through quantum gates, mathematically described as operators. 
			
			Take for example the Pauli operators. They can be written in terms of exterior products 
			using the computational basis as follows,
			
			\[
				\sigma_{0} = |0\rangle\langle0| + |1\rangle\langle1|
			\]
			\[
				\sigma_{x} = |0\rangle\langle1| + |1\rangle\langle0|
			\]
			\[
				\sigma_{y} = -i|0\rangle\langle1| + i|1\rangle\langle0|
			\]
			\[
				\sigma_{z} = |0\rangle\langle0| - |1\rangle\langle1|
			\]
			\\
			It is possible to write the above in OpenKet \footnote{This is not necessary 
			(though it helps a lot to explain how to build operators).
			The Pauli operators are defined inside OpenKet, as \texttt{X(\textsl{tag})}, 
			\texttt{Y(\textsl{tag})}, \texttt{Z(\textsl{tag})}, 
			($\sigma_{x}$, $\sigma_{y}$, $\sigma_{z}$, respectively). Where the \texttt{\textsl{tag}} 
			is an optional argument used to specify the operator tag, if handling more than 
			one qubit.}, having the advantage that they look visually the same,
			
			\texttt{
			\begin{tabbing}
				>>> sigma0 = Ket(0)*Bra(0) + Ket(1)*Bra(1); sigma0 \\
				|0><0| + |1><1| \\
				>>> sigmaX = Ket(0)*Bra(1) + Ket(1)*Bra(0); sigmaX \\
				|0><1| + |1><0| \\
				>>> sigmaY = -I*Ket(0)*Bra(1) + I*Ket(1)*Bra(0); sigmaY \\
				- I|0><1| + I|1><0| \\
				>>> sigmaZ = Ket(0)*Bra(0) - Ket(1)*Bra(1); sigmaZ \\
				|0><0| - |1><1| \\
			\end{tabbing}}
			
			They act on vectors,
			
			\texttt{
			\begin{tabbing}
				>>> v = 3*I*Ket(0) + 1.23*Ket(1); v \\
				3*I|0> + 1.23|1> \\
				>>> sigmaX*v \\
				1.23|0> + 3*I|1> \\
				>>> sigmaY*v \\
				- 1.23*I|0> - 3|1> \\ 
 				>>> sigmaZ*v \\
				3*I|0> - 1.23|1> \\
			\end{tabbing}}

			We can prove that $[\sigma_{x}, \sigma_{y}]=2i\sigma_{z}$,
			
			\texttt{
			\begin{tabbing}
				>>> left = sigmaX*sigmaY; left \\
				I|0><0| - I|1><1 \\
				>>> right = sigmaY*sigmaX; right \\
				- I|0><0| + I|1><1| \\
				>>> left - right \\
				2*I|0><0| - 2*I|1><1| \\		
			\end{tabbing}}
			
			And virtually any gate that we can represent as a sum of exterior products can be written 
			fairly easy using OpenKet.
		
\section{Built-in Functions}
	This functions are subroutines within OpenKet, so once you run the library all of them will be in 
	your workspace. 
	
	\subsection{Hermitian Conjugate: \texttt{Adj}}
		Gives the h.c. of the argument you input, be it a number (complex conjugate), a 
		\texttt{Ket}(\texttt{Bra}) (\texttt{Bra}(\texttt{Ket})), an \texttt{Operator} 
		(\texttt{AdjointOperator}) or a combination of them. For example,
		
		\texttt{
		\begin{tabbing}
			>>> Adj(I) \\
			-I \\
			>>> Adj(Ket(0)); Adj(Bra(1)) \\
			<0| \\
			|1> \\
			>>> Adj(I*Ket(0) - 4*Ket(1)) \\
			- I<0| - 4<1| \\
		\end{tabbing}}
			
		This is particularly useful for computing the density operator of a state, or its norm,
		
		\texttt{
		\begin{tabbing}
			>>> w = I*Ket(0) - 4*Ket(1); w \\
			I|0> - 4|1> \\
			>>> norm = Adj(w)*w; norm \\
			17 \\
			>>> rho = w*Adj(w); rho/norm \\
			1/17|0><0| - 4*I/17|0><1| + 4*I/17|1><0| + 16/17|1><1| \\
		\end{tabbing}}
		
	\subsection{Trace and Partial Trace: \texttt{Trace}, \texttt{TraceOut}}
		OpenKet allows us to trace or partially trace states expressed as a sum of exterior products. 
		Let's 
		first look into \texttt{TraceOut}. This has 2 arguments: the first one is the total expression 
		and the second one is the operator tag of the space you will be tracing out.
		
		Let's first construct a vector, (2 qubit system, part A and part B) and its density operator,
		
		\texttt{
		\begin{tabbing}
			>>> P = Ket(0, "A")*Ket(1, "B") - Ket(1, "A")*Ket(0, "B"); P \\
			|0\_A>|1\_B> - |1\_A>|0\_B> \\
			>>> R = P*Adj(P); R \\
			|0\_A>|1\_B><1\_B|<0\_A| - |0\_A>|1\_B><0\_B|<1\_A| - \\
			|1\_A>|0\_B><1\_B|<0\_A| + |1\_A>|0\_B><0\_B|<1\_A| \\
		\end{tabbing}}  
		
		So tracing out looks like this,
		
		\texttt{
		\begin{tabbing}
			>>> RB = TraceOut(R, "A"); RB \\
			|0\_B><0\_B| + |1\_B><1\_B| \\
			>>> RA = TraceOut(R, "B"); RA \\
			|0\_A><0\_A| + |1\_A><1\_A| \\
		\end{tabbing}}
		
		The total trace is done with the function \texttt{Trace}\footnote{A second optional argument 
		is the base which you would like to trace with, given as a list} . You only need to input the 
		total expression, \texttt{Trace} first finds all the operator tags and will then partially trace 
		them all out. 
		
		\texttt{
		\begin{tabbing}
			>>> Trace(R) \\
			2 \\
			>>> Trace(Ket(0)*Bra(1) + Ket(0)*Bra(0) + 2*I*Ket(1)*Bra(1)) \\
			1 + 2*I \\
		\end{tabbing}}
		
		One application of this is to be able to find if we are dealing with a mixed or a pure state,
		by computing the square of the density operator and calculating its total trace.
		
		Suppose we have the pure state \texttt{R1}= $|\psi\rangle\langle\psi|$ with, 
		$|\psi\rangle = (|0\rangle - |1\rangle)/\sqrt{2}$ and we have 
		the state \texttt{R2}= $0.25|0\rangle\langle0| + 0.75|1\rangle\langle1|$. First we 
		input this states,
		
		\texttt{
		\begin{tabbing}
			>>> psi = (Ket(0) - Ket(1))/sqrt(2) \\
			>>> R1 = psi*Adj(psi); R1 \\
			0.5|0><0| - 0.5|0><1| - 0.5|1><0| + 0.5|1><1| \\
			>>> R2 = 0.25*Ket(0)*Bra(0) + 0.75*Ket(1)*Bra(1); R2 \\
			0.25|0><0| + 0.75|1><1| \\
		\end{tabbing}}
		
		Both of them have total traces equal to one, but if we square them and calculate their traces,
		
		\texttt{
		\begin{tabbing}
			>>> Trace(R1*R1) \\
			1.0 \\
			>>> Trace(R2*R2) \\
			0.625 \\
		\end{tabbing}}
		
		Thus concluding that \texttt{R2} is not a pure state.
		
	\subsection{Matrix Representation: \texttt{Qmatrix}}
		We are able nonetheless to work with the matrix representation of operators. This is done 
		via \texttt{Qmatrix}. The function has 2 arguments, the first is the expression (sum of 
		exterior products) and the second is the base in list form. For example, recall the Pauli 
		operators of section \ref{quantumgates}; let's obtain the matrix representation of 
		$\sigma_{x}$ in the computational basis $\{|0\rangle, |1\rangle\}$,
		
		\texttt{
		\begin{tabbing}
			>>> basecomp = [Ket(0), Ket(1)] \\
			>>> Qmatrix(sigmaX, base) \\
			\lbrack0, 1] \\
			\lbrack1, 0] \\
		\end{tabbing}}
		
		Yet if we change base, for example $\{|+\rangle, |-\rangle\}$, where 
		$|\pm\rangle = \left( |0\rangle \pm |1\rangle \right) / \sqrt{2}$, the matrix representation 
		changes accordingly,
		
		\texttt{
		\begin{tabbing}
			>>> basepm = [(Ket(0)+Ket(1))/sqrt(2), (Ket(0)-Ket(1))/sqrt(2)] \\
			>>> Qmatrix(sigmaX, basepm) \\
			\lbrack1.0,\qquad0] \\
			\lbrack\quad0, -1.0] \\
		\end{tabbing}}
		
		This is useful when calculating for example eigenvalues, making things easier when computing 
		more complex properties such as Von Neumann entropy.		
	
\section{Time Evolution Examples}
	\subsection{The Problem}
		One important aspect of understanding a certain system is to understand its time evolution. 
		In quantum mechanics there are many ways of representing how a system evolves. In 
		this section we will demonstrate the use of OpenKet subroutines to solve numerically this 
		problem by considering first the Von Neumann equation,
		
		\[
			\dot{\rho} = -\frac{i}{\hbar}[\hat{H},\rho]
		\]
		\\
		Numerically this is a set of coupled ODE's (ordinary differential equations). 
		The thing is, one has to find a way to write 
		and name the variables to produce a script readable by an ODE's solver. 
		
		The problems presented in this manual are discrete level atoms in special situations. 
		We face this type of problems by first writing $\hat{H}$ as a sum 
		of exterior products of elements of a base; then verifying the time evolution of the matrix 
		elements of $\rho$, namely $\langle i|\rho|j\rangle$, considering them as variables.
		
		It is convenient to end up with real coefficients equations, due to the fact that not all 
		ODE's solvers work with complex numbers, and the point is to make this as transparent as possible.
	
	\subsection{The Translator \texttt{QeqN}}
		The first thing is to write the Hamiltonian and every extra term in the right hand side of 
		the time evolution equation in terms that OpenKet understands. Then comes in the function 
		\texttt{QeqN}. Basically it translates terms of the form $\langle i|\rho|j\rangle$ to 
		symbolic expressions (variables) which are easier to handle. Then it writes the equations into 
		a text file as well as the dictionary relating the $\langle i|\rho|j\rangle$ terms with the 
		variables. 
		
		\texttt{QeqN} has the following arguments
		\texttt{QeqN}\textsl{(density operator, time derivative, base, file name, dictionary name)}. 
		The fist argument is just an object class \texttt{Operator} used to designate the density 
		operator. The second is the right hand side of the time evolution equation. The third argument 
		is a list, containing the elements of the base, expressed as kets. The fourth and fifth 
		arguments are strings naming the text file and the dictionary written in the file. 
		
		\texttt{QeqN} uses 2 other functions: \texttt{Dictionary}\textsl{(base, density operator)}
		\footnote{Not to be confused with the dictionary written in the file} 
		and \texttt{Qch}\textsl{(expression, dictionary)}. 
		
		\texttt{Dictionary} does just that,
		create a dictionary (with name D) with words as $\langle i|\rho|j\rangle$ objects and 
		definitions as variables. It is worth noting that in order to have always real 
		coefficients, the variables must be in fact 
		$\frac{\langle i|\rho|j\rangle + \langle j|\rho|i\rangle}{2}$ and 
		$\frac{\langle i|\rho|j\rangle - \langle j|\rho|i\rangle}{2i}$; 
		because $\langle i|\rho|j\rangle$ is 
		in general complex\footnote{Although $\langle i|\rho|i\rangle$ is always real in physical 
		situations}. Let's illustrate this by building the dictionary of the base 
		$\{|0\rangle, |1\rangle\}$, 
		
		\texttt{
		\begin{tabbing}
			>>> base=[Ket(0), Ket(1)]; R=Operator("R") \\
			>>> Dictionary(base, R) \\
			\{ '<0|R|0>': Re0, \\
 			\quad'<0|R|1>': Re1 + I*Im1, \\
 			\quad'<1|R|0>': Re1 - I*Im1, \\
 			\quad'<1|R|1>': Re2\} \\
		\end{tabbing}}
		
		The symbols \texttt{Re0, Re1, Im1, Re2} are declared as real variables.
		
		The function \texttt{Qch} uses the dictionary created above to substitute bracket terms, for 
		the symbols declared in \texttt{Dictionary}. For example,
		
		\texttt{
		\begin{tabbing}
			>>> D = Dictionary(base, R) \\
			>>> expr = Bra(0)*R*Ket(0)+2*I*Bra(0)*R*Ket(1)-3.4*Bra(1)*R*Ket(1); expr \\
			<0|R|0> + 2*I<0|R|1> - 3.4<1|R|1> \\
			>>> EXPR = Qch(expr, D); EXPR \\
			Re0 - 3.4*Re2 + 2*I*(Re1 + I*Im1) \\
		\end{tabbing}}
		
		So \texttt{QeqN} first creates a dictionary, then substitutes bracket terms, and then writes to 
		a file. Before writing to the file, it replaces the real variables with strings of the form 
		\texttt{y[n]} where n is some number. The cost of introducing new notation is justified because 
		strings are objects which consume a lot less memory than variables, in this case.\footnote{Why 
		then declare them first as real variables? Because we need to be able to operate 
		and simplify expressions, and this is done on symbolic variables, not strings. 
		The package Sympy takes care of all the algebra, but again, in order to do so, the objects 
		must be symbolic variables.}
		The written file looks something like this,
		
		\texttt{
		\begin{tabbing}
			def f(y, t): \\
    			\qquad return [-2*y[2],0,y[0] - 1.0*y[3],2*y[2]] \\
			dict=\{'<1|R|0>': y1 - I*y2, '<1|R|1>': y3, \\
			'<0|R|0>': y0, '<0|R|1>': y1 + I*y2\} \\
		\end{tabbing}}
		
		The text is the usual way to define a function in Python. The name \texttt{f} is always used 
		as well as its arguments. 
		
		\texttt{dict} is a dictionary object, having its name specified in the last argument of 
		\texttt{QeqN}. It serves only the purpose of aiding the user in determining which bracket term 
		is which string\footnote{Helping in the study of observables, as will be evident in the last 
		section}.
		
		After the \texttt{return} command, there is a list which contains the temporal evolution of 
		the matrix elements. The position of a certain \texttt{item} determines the time derivative of 
		the matrix element corresponding to \texttt{y[item]}. For example, in the text presented before, 
		we have \texttt{[-2*y[2],0,y[0] - 1.0*y[3],2*y[2]]}, so looking at the first entry,
		\\
		
		$\frac{d}{dt}\texttt{y[0]} = -2\texttt{y[2]}$
		\\
		
		Equivalent to (using the \texttt{dict}),
		
		\[
			\frac{d}{dt}\langle0|\rho|0\rangle = 
			-2\frac{\langle0|\rho|1\rangle - \langle1|\rho|0\rangle}{2i}
		\]
		
		The sintaxis of the file is such that the ODE solver \texttt{odeint} from the package Scipy 
		understands it. So all we need to do, is run the file from our workspace, and specify the 
		time step and the initial condition. To illustrate this, we show two interesting
		physical problems.
					
	\subsection{Two Level Semiclassical Atoms}
		In some physical situations when an atom interacts with an electromagnetic field, we can 
		describe the atom using a 2 element base; let this base be $\{|0\rangle, |1\rangle\}$. 
		Furthermore let the energy of the levels be $\hbar\omega_{0}$ and $\hbar\omega_{1}$ 
		respectively. Now, the total Hamiltonian may be divided into two parts, the Hamiltonian 
		of the unperturbed atom $\hat{H}_{0}$, and the interaction part $\hat{V}$. These can be 
		found to have the following form,
		
		\[
			\hat{H}_{0} = \hbar\omega_{0}|0\rangle\langle0| + 
			\hbar\omega_{1}|1\rangle\langle1| \\
		\]
		\[
			\hat{V} = 
			gE(t)|0\rangle\langle1| + 
			g^{\ast}E(t)|1\rangle\langle0|
		\]
		\\
		Where we have put $g = e\langle0|\vec{r}|1\rangle \cdot \vec{\epsilon}$ and 
		$\vec{E}(t) = |\vec{E}(t)|\vec{\epsilon} \equiv E(t)\vec{\epsilon}$ with 
		$|\vec\epsilon|=1$. We assume the field has the form,
		
		\[
			E(t) = A\cos(\nu t)
		\]
		
		It is convenient to change to the interaction picture, and following the 
		rotating wave approximation (RWA) we end up with the Hamiltonian (setting 
		$\Omega = \frac{gA}{2\hbar}$ and having the frequency of the field equal to the difference 
		of the two level frequencies),
		
		\[
			\hat{H} = 
			\hbar( \Omega|0\rangle\langle1| +
			\Omega^{\ast}|0\rangle\langle1|)
		\]
		
		Following is the script which numerically solves the time evolution of 
		the system described\footnote{assuming $\Omega=1$ for simplicity},
		
		\texttt{
		\begin{tabbing}
			run -i openket.py \\
			from scipy.integrate import odeint 
			\textsl{\# import Scipy's ODE solver}\\
			R = Operator("R") \\
			b = [Ket(0), Ket(1)] \\
			H =  Ket(0)*Bra(1) + Ket(1)*Bra(0) \\
			Rdot = -I*Commutator(H, R) \\
			QeqN(R, Rdot, b, "func", "dic") \\
			run -i func \\
			initialc = [1,0,0,0] 
			\textsl{\# initial condition}\\
			t = linspace(0,10,1000) 
			\textsl{\# time-step}\\
			sol = odeint(f, initialc, t) \\
			es = sol[:,3] \\
			plot(t, es) \\
		\end{tabbing}}
		
		The graph we obtain is the following,
		\newpage
		\begin{figure}[h]
			\centering\includegraphics[scale=.40]{manual1.pdf}
			\caption{Time evolution for $\langle1|\rho|1\rangle$ with the initial condition 
			$\rho_{(t=0)} = |0\rangle\langle0|$}
		\end{figure}
		
		Note that the initial condition is given as a list, where the position in the list indicates 
		the bracket term that one wishes to give as the initial condition (the \texttt{dic} tells 
		you which is which). Also the time-step is made with the command \texttt{linspace} having 
		3 arguments: the first two specify the interval of solution while the third one is the total 
		number of points in the interval. In the script, \texttt{sol} is the actual solution, and has the 
		form,
		
		\[
			\texttt{sol} = [\overbrace{[\texttt{y[0]}, \texttt{y[1]}, \texttt{y[2]}, 
			\texttt{y[3]} ]}^{t=0}, \overbrace{[\texttt{y[0]}, \texttt{y[1]}, \texttt{y[2]}, 
			\texttt{y[3]} ]}^{t=0.01}, \dots, \overbrace{[\texttt{y[0]}, \texttt{y[1]}, 
			\texttt{y[2]}, \texttt{y[3]} ]}^{t=10}]
		\]
		
		So the command \texttt{es = sol[:,3]} partitions the total solution \texttt{sol} and only 
		keeps the third entry (\texttt{y[3]} or $\langle1|\rho|1\rangle$) in all the lists within 
		\texttt{sol}. The last line is just the plotting of the probability of finding the atom in the
		excited state versus time.
		\\
		
		Suppose we don't want to make the RWA, or treat the problem from the interaction picture; 
		that is, work with the total Hamiltonian\footnote{Note that it is time-dependent explicitly}
		$\hat{H}=\hat{H}_{0}+\hat{V}$,
		
		\[
			\hat{H} = \hbar\omega_{0}|0\rangle\langle0| + 
			\hbar\omega_{1}|1\rangle\langle1| + 
			2\hbar\Omega\cos(\nu t)|0\rangle\langle1| + 
			2\hbar\Omega^{\ast}\cos(\nu t)|1\rangle\langle0|
		\]
		
		A possible script, (notice that the structure is the same as the last script)\footnote{Setting 
		$\Omega = 0.5$, again for simplicity}:
		\newpage
		
		\texttt{
		\begin{tabbing}
			run -i openket.py \\
			from scipy.integrate import odeint \\
			from sympy import * 
			\textsl{\# the whole package} \\
			\textsl{\# to be able to treat cos(n*t) as a variable} \\
			R = Operator("R") \\
			b = [Ket(0), Ket(1)] \\
			w1=10; w0=1; v=9 
			\textsl{\# define parameters (resonant condition: v = w1-w0)} \\
			t = var("t") \\
			H = w0*Ket(0)*Bra(0) + w1*Ket(1)*Bra(1) \\
			H = H + cos(v*t)*(Ket(0)*Bra(1) + Ket(1)*Bra(0)) \\
			Rdot = -I*Commutator(H, R) \\
			QeqN(R, Rdot, b, "func", "dic") \\
			run -i func \\
			initialc = [1,0,0,0] \\
			t = linspace(0,10,1000) 
			\textsl{\# has to have the same tag as the variable t}\\
			sol = odeint(f, initialc, t) \\
			es = sol[:,3] \\
			plot(t, es) \\
		\end{tabbing}}
		
		We obtain,
		
		\begin{figure}[h]
			\centering\includegraphics[scale=.40]{manual2.pdf}
			\caption{Time evolution for $\langle1|\rho|1\rangle$ with the initial condition 
			$\rho_{(t=0)} = |0\rangle\langle0|$ and without the RWA. Values 
			($\omega_{0}=1; \omega_{1}=10; \nu=9$)}
		\end{figure}
		
		Changing the parameters a bit,
		
		\begin{figure}[h]
  			\centering
  			\subfloat[Full size]{\includegraphics[scale=.35]{manual3.pdf}}                
 			\subfloat[Detail]{\includegraphics[scale=.35]{manual4.pdf}}
  			\caption{Time evolution for $\langle1|\rho|1\rangle$ with the initial condition 
			$\rho_{(t=0)} = |0\rangle\langle0|$ and without the RWA. Values 
			($\omega_{0}=1; \omega_{1}=100; \nu=99$)}
		\end{figure}
				
	\subsection{Two Level Atom Inside a Cavity}
		The last scenario tests basically every aspect of OpenKet. Let a two level atom inside a 
		cavity interact with an electromagnetic field. There is laser pumping and losses, and there is 
		spontaneous emission from the atom. This can be modeled by having a Hamiltonian of the form,
		
		\[
			\hat{H} = 
			g\left( \hat{a}^{\dagger}\sigma^{-} + 
			\hat{a}\sigma^{+} \right) + 
			\epsilon \left( \hat{a}^{\dagger} + 
			\hat{a} \right)
		\]
	
		Where $\sigma^{-}=|0\rangle\langle1|; \sigma^{+}=|1\rangle\langle0|$. The Von Neumann equation 
		has extra dissipative terms,
		
		\[
			\dot{\rho} = 
			- \frac{i}{\hbar}\left[\hat{H}, \rho\right] +
			kL(\hat{a}) + 
			\frac{\gamma}{2}L(\sigma^{-})
		\]

		Having $L(\hat{o}) = 2\hat{o}\rho\hat{o}^{\dagger} - \rho\hat{o}^{\dagger}\hat{o} - 
		\hat{o}^{\dagger}\hat{o}\rho$. 

		To solve this we have to use a finite base (truncating the Fock states), such as,
		
		\[
			\{|i\rangle_{atom}|j\rangle_{photon} : i=0,1\quad;\quad j=0,1,\dots,n\}
		\]

		Deciding which $n$ has to be carefully thought, (both because of physical issues, and 
		computer speed issues). 
		
		Let us work with a base having this structure and with $n=4$. Writing the script is no problem, 
		it's just copying the Hamiltonian and time evolution equations from before
		\footnote{the base will have 10 elements 
		$\{|0\rangle_{at}|0\rangle_{ph}, |0\rangle_{at}|1\rangle_{ph}, \dots,
		|1\rangle_{at}|0\rangle_{ph}, |1\rangle_{at}|1\rangle_{ph}, \dots, 
		|1\rangle_{at}|4\rangle_{ph}\}$ producing a 100 element density matrix. When we input the 
		initial condition, it must be in the form of a list of 100 entries. This is not difficult at 
		all, as may be evident from the script. The challenge is to obtain a certain observable.}
		\newpage
		
		\texttt{
		\begin{tabbing}
			run -i openket.py \\
			from scipy.integrate import odeint \\
			a = AnnihilationOperator("photon",4) 
			\textsl{\# maximum value 4} \\
			aa = CreationOperator("photon",4)
			\textsl{\# maximum value 4} \\
			splus = Ket(1,"atom")*Bra(0,"atom") \\
			sminus = Adj(splus) \\
			g = 0.6; epsilon = 0.1; k = 0.5; gamma = 1.5
			\textsl{\# define parameters} \\
			H = g*(aa*sminus + a*splus) + epsilon*(a + aa)
			\textsl{\# Hamiltonian} \\
			def L(x): \textsl{\# Define the function L} \\
			\qquad	R=Operator("R") \\
			\qquad	return 2*x*R*Adj(x) - R*Adj(x)*x - Adj(x)*x*R \\
			R = Operator("R") \\
			Rdot = - I*Commutator(H, R) + \textbackslash \\
			k*(2*a*R*aa - R*aa*a - aa*a*R) + (gamma/2.0)*L(sminus) \\
			l=[]; m=[]; b=[]  
			\textsl{\# start the base construction} \\
			for i in range(2): \\
			\qquad    l.append(Ket(i, "atom")) \textsl{\# atom part} \\
			for j in range(5): \\ 
			\qquad    m.append(Ket(j,"photon")) \textsl{\# photon part} \\
			for i in l: \\
 			\qquad   for j in m: \\
			\qquad \qquad	    b.append(i*j) \textsl{\# put them together} \\
			y0 = [] \textsl{\# initial condition} \\
			for i in range(100): \\
			\qquad	y0.append(0) \textsl{\# fill it with zeroes} \\
			y0[0] = 1 \textsl{\# at time t=0 we have atom at base state and no photons} \\
			t = linspace(0,10,100) \\
			QeqN(R, Rdot, b, "func", "dic") \\
			run -i func \\
			sol = odeint(f, y0, t) \\
		\end{tabbing}}

		This completely solves the problem numerically. \texttt{sol} has all the time evolution 
		information of $\rho$. The problem is, how can we manage this information in order to say, 
		plot the time evolution of a certain observable?
		
		\subsubsection{Observables and \texttt{SubsSol}}	

			Suppose now we want to plot the mean number of photons 
			against time. The observable is $tr(\rho\hat{a}^{\dagger}\hat{a})$.
			To be able to plot this, we must first compute the trace and have something of the type,
			$\sum_{i,j}c_{ij}\langle i|\rho|j\rangle$ This can be done directly with the function 
			\texttt{Trace} using the second optional argument, (keeping the names we defined in 
			the previous script)
		
			\texttt{
			\begin{tabbing}
				>>> obs = Trace(R*aa*a, b); obs \\
				<1\_photon|<0\_atom|R|0\_atom>|1\_photon> + <1\_photon|<1\_atom|R|1\_atom>|1\_photon> + \\
				2<2\_photon|<0\_atom|R|0\_atom>|2\_photon> + 2<2\_photon|<1\_atom|R|1\_atom>|2\_photon> + \\
				3<3\_photon|<0\_atom|R|0\_atom>|3\_photon> + 3<3\_photon|<1\_atom|R|1\_atom>|3\_photon> + \\
				4<4\_photon|<0\_atom|R|0\_atom>|4\_photon> + 4<4\_photon|<1\_atom|R|1\_atom>|4\_photon>\\
			\end{tabbing}}
		
			The point is to transform the above expression into a symbolic variable type. Remember the 
			dictionary object \texttt{dic}, (obtained from \texttt{QeqN})
		
			\texttt{
			\begin{tabbing}
				>>> dic \\
				\{ '<0\_photon|<0\_atom|R|0\_atom>|0\_photon>': y0, \\
 				\quad'<0\_photon|<0\_atom|R|0\_atom>|1\_photon>': y1 + I*y2,\\
 				\quad'<0\_photon|<0\_atom|R|0\_atom>|2\_photon>': y3 + I*y4,\\
				\quad\dots \\
				\quad'<0\_photon|<1\_atom|R|0\_atom>|2\_photon>': y41 - I*y42, \\
 				\quad'<0\_photon|<1\_atom|R|0\_atom>|3\_photon>': y54 - I*y55, \\
 				\quad'<0\_photon|<1\_atom|R|0\_atom>|4\_photon>': y65 - I*y66, \\
 				\quad'<0\_photon|<1\_atom|R|1\_atom>|0\_photon>': y75, \\
 				\quad'<0\_photon|<1\_atom|R|1\_atom>|1\_photon>': y76 + I*y77, \\
				\quad\dots \} \\
			\end{tabbing}}

			So using this dictionary and the function \texttt{Qch}, we are able to 
			transform the expression \texttt{obs},
		
			\texttt{
			\begin{tabbing}
				>>> obs = Qch(obs, dic); obs \\
				y19 + y84 + 2*y36 + 2*y91 + 3*y51 + 3*y96 + 4*y64 + 4*y99\\
			\end{tabbing}}
			
			Now comes \texttt{SubsSol}. At this point, we have the instructions of what we want to plot (the 
			above expression) and the whole time evolution (the object \texttt{sol}); to plot this manually 
			would be very unpractical (you would have to type \texttt{sol[:,19] + sol[:,84] + 
			2*sol[:,36] +} \dots). 
			
			\texttt{SubsSol} does this for you\footnote{it also supports expressions
			with functions in them e.g. \textsl{sin, log, cos}, etc. (actually all functions
			recognized by Sympy)}. It has two arguments 
			\texttt{SubsSol}(\textsl{obs, sol}). The first argument is the instructions of what to plot, 
			(the expression with the \texttt{yi's} variables); and the second argument is the object \texttt{sol}
			which returns \texttt{odeint}, i.e. the whole time evolution. It looks like this,
			
			\texttt{
			\begin{tabbing}
				>>> OBS = SubsSol(obs, sol)
			\end{tabbing}}
			
			Then \texttt{OBS} is an array of elements which has the time evolution of the mean number of photons,
			in this case. Plotting \texttt{OBS} results in,
			
			\begin{figure}[h]
				\centering\includegraphics[scale=.40]{manual5.pdf}
				\caption{Evolution of the mean number of photons, having the atom at the base state
				and no photons as an initial condition}
			\end{figure}

			Let's try now plotting the time evolution of the probability of finding the atom at the base state.
			This can be calculated with $tr(\rho|0\rangle\langle0|_{atom}\otimes\mathbb{1}_{photon})$. Typing this in 
			OpenKet is easy:
			
			\texttt{
			\begin{tabbing}
				>>> bs=Trace(R*Ket(0, 'atom')*Bra(0, 'atom'), b); bs \\
				<0\_photon|<0\_atom|R|0\_atom>|0\_photon> + <1\_photon|<0\_atom|R|0\_atom>|1\_photon> + \\
				<2\_photon|<0\_atom|R|0\_atom>|2\_photon> + <3\_photon|<0\_atom|R|0\_atom>|3\_photon> + \\
				<4\_photon|<0\_atom|R|0\_atom>|4\_photon> \\
			\end{tabbing}}

			Note that the tensor product with the identity on the photon part is implicit. Next, we change it to 
			sympy variables using the \texttt{dic},
			
			\texttt{
			\begin{tabbing}
				>>> bs = Qch(bs, dic); bs \\
				y0 + y19 + y36 + y51 + y64\\
			\end{tabbing}}

			Obtain the numerical information with \texttt{SubsSol},
			
			\texttt{
			\begin{tabbing}
				>>> BS = SubsSol(bs, sol)\\
				>>> plot(t, BS)
			\end{tabbing}}

			The plot\footnote{This can be done in one line, due to the fact that Python allows for nesting of 
			functions: \newline \texttt{>>> plot(t, SubsSol(Qch(Trace(R*Ket(0,'atom')*Bra(0,'atom'), b), dic), sol, 
			b))}},
			
			\begin{figure}[h]
				\centering\includegraphics[scale=.40]{manual6.pdf}
				\caption{Time-evolution of the Probability of finding the atom in the base state; 
				having the atom at the base state and no photons as an initial condition}
			\end{figure}
			
			


















\end{document}